Defining a custom bounding box

Displaying (and some computations) the entire network we just downloaded can slow down R environment. For class exercise, we will do the calculation for the whole network but limit the display to a smaller area. For that, we will create a custom bbox using coordinates of our selection. This technique of generating your own bounding box can also be useful if you have a specific area of interest that doesn’t overlap well with commonly used boundaries.

You need to define two points for bounding box, one point at the lower left corner and the other at the upper right corner. You can go to Google Maps, right-click on a point on map, and copy the XY coordinate. Let’s store them in p1 and p2.

# p1 is lower left corner, p2 is the upper right corner
p1 <- c(33.746217847959734, -84.40851957882589)
p2 <- c(33.785889694219634, -84.36354430149285)

You will then need to format the two coordinates into the same format as bb, the bounding box object we created above. Then, let’s convert it into sf object.

# Custom BB
my_bb <- matrix(c(p1[2], p1[1],
                  p2[2], p2[1]), ncol = 2)

colnames(my_bb) <- c("min", "max")
rownames(my_bb) <- c("x", "y")

# Custom BB to sf
my_bb_sf <- my_bb %>% t %>% data.frame() %>% 
  st_as_sf(coords = c("x", "y"), crs = 4326) %>% 
  st_bbox() %>% 
  st_as_sfc()

## Plot--
tmap_mode('view')

## tmap mode set to interactive viewing

tm_shape(bb_sf) + tm_borders(col = "black") +  # Black = larger bbox
  tm_shape(my_bb_sf) + tm_borders(col = "red")  # Red = smaller bbox

Simplifying network

Multiple edges can connect the same pair of nodes, called multiple edges. There can also be loops that starts and ends at the same node (e.g., think of a cul-de-sac). The former case can be detected using edge_is_multiple() and the latter using edge_is_loop().

When removing a set of multiple edges using functions shown below, they always keep the first edge and discard others. By arranging the order of edges for each set of multiple edges, you can specify which edge you want to preserve.

This way of simplification means that, when the multiple edges within a set have different attributes, all attribute information except for the preserved one would be lost. In such cases, you can merge those edges. Then, the output would have the geometry of the first edge in the set, but the attributes would be some summary of all the edges in the set. to_spatial_simple() function does this work.

# Let's simplify our network
simple_net <- net %>%
  activate("edges") %>%
  filter(!edge_is_multiple()) %>%
  filter(!edge_is_loop()) 

# Because the difference is not really discernible, we just print out the differences in the edge count.
message(str_c("Before simplification, there were ", net %>% st_as_sf("edges") %>% nrow(), " edges. \n",
            "After simplification, there are ", simple_net %>% st_as_sf("edges") %>% nrow(), " edges."))

## Before simplification, there were 1609 edges. 
## After simplification, there are 1600 edges.

Subdivide edges

When as_sfnetwork() function converts an sf linestrings, the nodes are defined as the endpoints of each linestring. If you zoom into Midtown in the map above, you will see that there are many intersections that do not have nodes, which are errors. We can use to_spatial_subdivision predicate to create points at the intersections and cut the edges accordingly.

To better understand

# Subdivision
subdiv_net <- convert(simple_net, sfnetworks::to_spatial_subdivision)

## Warning: to_spatial_subdivision assumes attributes are constant over geometries

# Add a custom index -> This will be used for visualization purposes later
subdiv_net <- subdiv_net %>% 
  activate("nodes") %>% 
  mutate(custom_id = seq(1, subdiv_net %>% st_as_sf("nodes") %>% nrow()),
         is.new = case_when(is.na(.tidygraph_node_index) ~ "new nodes",
                            TRUE ~ "existing nodes"),
         is.new = factor(is.new)) 

## Plot--
subdiv_pal <- colorFactor(palette = c("yellow", "red"), domain = subdiv_net %>% st_as_sf("nodes") %>% pull(is.new))

subdiv_map <- leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  setView( -84.3854, 33.7668, zoom = 14) %>% # zooming in to show more details
  addPolylines(data = subdiv_net %>% st_as_sf('edges'), col = "grey", weight = 3, opacity = 0.9) %>% 
  addCircles(data = subdiv_net %>% st_as_sf('nodes'), fillColor = ~subdiv_pal(is.new), stroke = F, radius = 20, fillOpacity = 0.7) %>% 
  addLegend(position = "bottomright", pal = subdiv_pal, values = subdiv_net %>% st_as_sf("nodes") %>% pull(is.new)) 

subdiv_map

Delete pseudo nodes

smoothed_net <- convert(subdiv_net, sfnetworks::to_spatial_smooth)

# Extract removed points
removed <- setdiff(subdiv_net %>% st_as_sf('nodes') %>% pull(custom_id), 
                   smoothed_net %>% st_as_sf('nodes') %>% pull(custom_id))

smooth_map <- leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  setView( -84.3854, 33.7668, zoom = 14) %>% # zooming in to show more details
  addPolylines(data = smoothed_net %>% st_as_sf('edges'), col = "grey", weight = 3, opacity = 0.9) %>% 
    addCircles(data = smoothed_net %>% st_as_sf('nodes'), fillColor = 'yellow', stroke = F, radius = 20, fillOpacity = 0.7) %>% 
    addCircles(data = subdiv_net %>% st_as_sf("nodes") %>% filter(custom_id %in% removed), 
               fillColor = "red", stroke = F, radius = 15, fillOpacity = 0.8) %>% 
  addControl(html = htmltools::HTML("<b>Red points denote removed nodes</b>"), position = "bottomright")

smooth_map

Shortest path

# Change crs for convenience
smoothed_net <- smoothed_net %>% st_transform(4326)

# Start point
start_p <- st_point(c(-84.40364459476174,33.776160322717544)) %>% st_sfc(crs = 4326) # CRC at Georgia Tech
# End point
target_p <- st_point(c(-84.37639335217811, 33.75718076235044)) %>% st_sfc(crs = 4326) # MLK National Historical Park
# Get the shortest path
paths = st_network_paths(smoothed_net, from = start_p, to = target_p)
# Extract shortest path
paths_sf <- smoothed_net %>%
  slice(paths$node_paths[[1]])

# Visualize
leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  setView( -84.3854, 33.7668, zoom = 14) %>% # zooming in to show more details
  addPolylines(data = smoothed_net %>% st_as_sf("edges"), color = 'grey', weight = 2, opacity = 0.7) %>% 
  addPolylines(data = paths_sf %>% st_as_sf("edges"), color = "red", weight = 4, opacity = 0.8) %>% 
  addCircles(data = paths_sf %>% st_as_sf("nodes"), stroke = F, fillColor = "red", fillOpacity = 0.8, radius = 50)

## Assignment: sample one point from each tract, calculate average travel time from one census tract to all others. Repeat for all tracts.
## See if there is any patterns to the accessibility.

Extract intersections

end_points <- smoothed_net %>% 
  st_as_sf('nodes') %>% 
  st_join(smoothed_net %>% activate("edges") %>% st_as_sf())

intersections <- end_points %>% 
  group_by(.tidygraph_node_index) %>% 
  summarise(n = n()) %>% 
  filter(n > 1) 

leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  setView( -84.3854, 33.7668, zoom = 14) %>% # zooming in to show more details
  addPolylines(data = smoothed_net %>% st_as_sf('edges'), col = "grey", weight = 3, opacity = 0.9) %>% 
    addCircles(data = intersections %>% st_as_sf('nodes'), fillColor = 'red', stroke = F, radius = 20, fillOpacity = 0.7) 

Spatial join the Census data with our network

# Before we do the join, let's use centrality measures that were calculated on the entire network.
# It takes time, so I calculated it and saved it.
# You can read the data back in
load("D:/Dropbox (GaTech)/Work/Working/School/UA_2022/external/Lab/module_2/network_char.RData")


census_centrality <- acs2020c %>%
  st_join(network_char %>% st_as_sf("nodes") %>% st_transform(crs = 4326)) %>% # When I calculated centrality measures, the CRS was 26967, so reverting it back to 4326
  group_by(GEOID) %>%
  summarise(n = n(),
            hhinc = mean(hhinc, na.rm = T),
            pct_wh = mean(pct_wh, na.rm = T),
            pct_bl = mean(pct_bl, na.rm = T),
            pct_novhc = mean(pct_novhc, na.rm = T),
            node_bc = mean(node_bc, na.rm = T))

tm_shape(acs2020c) + tm_polygons(col = "grey", alpha = 0.5) +
  tm_shape(census_centrality) + tm_polygons(col = "node_bc", style = "quantile") +
  tm_shape(bb_sf) + tm_borders()

census_centrality_plot <- census_centrality %>%
  mutate(hhinc = log(hhinc),
         pct_novhc = log(pct_novhc + 0.02)) %>%
  pivot_longer(cols = c('hhinc', 'pct_wh', 'pct_bl', 'pct_novhc'), names_to = "variable", values_to = "value") %>% 
  mutate(variable = factor(variable, labels = c("Household Income (logged)", "% Black", "% No HH with no cars (logged)", "% White")))

centrality_plot <- census_centrality_plot %>%
  ggplot() +
  geom_point(aes(x = node_bc, y = value), alpha = 0.2) +
  facet_wrap(~variable, scales = "free_y") +
  labs(x = "Centrality", title = "Centrality VS. Socio-demographics") +
  theme_bw()

centrality_plot + ggpubr::stat_cor(aes(x = node_bc, y = value))

## Warning: Removed 896 rows containing non-finite values (stat_cor).

## Warning: Removed 896 rows containing missing values (geom_point).

ggplotly(centrality_plot)